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SUMMARY 


In  reference  to  any  solution  of  a  conservative 
dynamical  system  with  two  degrees  of  freedom,  Hill's  equation 
is  generalized  to  encompass  non-necessarily  isoenergetic 
displacements  as  well  as  the  isoenergetic  displacements 
caused  by  a  variation  of  a  parameter. 

This  new  variational  equation  is  made  the  foundation 
of  a  methodical  procedure  for  continuing  numerically  natural 
families  of  periodic  orbits.  The  method  consists  of  two  steps- 
an  isoenergetic  corrector  and  a  tangential  predictor. 

Although  the  algorithm  makes  no  assumption  of  symmetry 
on  the  periodic  orbits  to  be  continued,  special  attention  is 
paid  to  the  symmetric  orbits,  but  only  to  show  how  in  these 
cases  the  method  can  be  simplified  substantially. 


The  subject  matter  of  this  paper  is  nothing  else  but  the 
elementary  problem  of  continuing  numerically  analytic  manifolds  of 
periodic  orbits  for  a  conservative  dynamical  system  with  two  degrees 
of  freedom.  Apparently  for  the  first  time,  a  methodical  answer  is 
given  here.  Our  procedure,  which  throughout  adheres  strictly 
to  the  analytical  foundations  of  the  question,  aims  at  keeping  as 
close  as  possible  to  the  basic  technique  of  Poincard,  namely  the  method 
of  analytic  continuation  based  on  Cauchy's  local  existence  theorem. 

In  the  classical  contributions  to  this  problem  (Darwin,  Stromgren, 

Lemaitre) ,  the  main  issue  is  immediately  obscured  by  the  accidental 
fact  that,  for  some  families,  the  periodicity  conditions  can  be  replaced 
by  symmetry  conditions,  and  therefore  the  problem  of  correcting  initial 
conditions  is  replaced  by  that  of  adjusting  boundary  conditions. 

Evenmore,  in  recent  years,  in  spite  of  the  many  capabilities  offered 
by  electronic  computers,  the  numerical  continuation  of  a  manifold  of 
periodic  orbits  tended  to  degenerate  into  disreputable  tricks  based  on 
optical  illusions  rather  than  on  analytical  certainties.  The  result  has 
been  an  overabundance  ol'  numerical  material  whose  subjective  interpretation 
leads  to  conclusions  at  variance  with  propositions  firmly  established  by 
analysis.  The  classical  Instance  of  such  unfortunate  accidents  is  the 
still  open  controversy  concerning  a  genealogy  of  periodic  orbits 
established  numerically  by  Darwin  (1897)  and  questioned  by  Poincard  (1899) 
on  analytical  grounds. 


Of  course,  one  can  only  expect  to  find  here  the  most  elementary 
part  of  the  very  extensive  theory  of  periodic  orbits.  A  natural  family 
of  periodic  orbits  is  defined  by  the  local  existence  of  power  series  in 
Painleve"'s  constant  of  integration  to  represent  the  manifold  (§1).  The 
fundamental  result  about  a  periodic  orbit  is  an  extension  to  non- 
necessarily  isoenergetic  normal  displacements  of  Hill's  equation  which 
is  valid  only  for  isoenergetic  variations  (§2).  Therefrom  is  derived 
the  vital  possibility  of  converging  to  a  periodic  orbit  without  leaving 
the  energy  manifold  on  which  it  lies  (§4)  as  well  a  tangential  predictor 
(55)  which  expresses  the  "solidarity"  between  the  phase  states  of  a 
natural  family  on  different  periodic  orbits  of  the  manifold  where  it  is 
deprived  of  singularities.  Symmetric  periodic  orbits  are  given  here  a 
special  treatment  only  because  it  is  possible  to  obtain  them  by  integration 
over  only  half  an  estimate  of  their  period  (§7). 

After  the  two  fundamental  steps  of  our  numerical  continuation  have 
been  derived,  the  other  sections  of  this  paper  concern  what  is  probably 
the  most  useful  concept  in  the  theory  of  conservative  systems,  namely  the 
isoenergetic  rate  of  variation  of  the  state  variables  with  respect  to  a 
parameter.  We  show  how  such  rates  can  be  computed  intrinsically  here  again 
from  an  extension  of  Hill's  equation  (§8).  It  enables  us  to  transpose 
our  procedure  of  numerical  continuation  to  the  cases  when  new  time  variables 
are  introduced  for  any  fixed  value  of  the  energy  constant  (S9). 


1.  NUMERICAL  CONTINUATION  OF  A  NATURAL  FAMILY 


Given  a  conservative  dynamical  system  with  two  degrees  of  freedom 


£u  2^gllql  +  2g12qlq2  +  g22q2  *  +  flqi  +  f2q2  +  U 


one  can  always  choose  an  isothermal  set  of  coordinates  (; i,y )  and 
redefine  the  independent  variable  so  that  the  equations  of  motion  take 
on  the  simple  form 


x  *  2Ay  +  W  ,  y  *  -2Ax  +  W 

J  x  J  y 

where  A  and  W  are  functions  of  the  coordinates  x  and  y 
(Birkhoff  1915). 

The  equations  (1)  admit  the  integral 

C  -  2W  -  (x2  +  y2) 


which  it  will  be  convenient  to  refer  to  as  the  energy  integral;  C  is 
an  arbitrary  constant  of  integration,  to  be  called  the  constant  of  energy 
or  also  Painlevi' 8  constant  of  integration  (Chazy  1953). 

Let  us  assume  that,  for  the  initial  conditions  (xo*^o*^o'^0^  at 
time  t  ■  0,  the  equations  (1)  possess  a  solution 
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which  is  periodic.  We  denote  by  Tq  the  period  of  this  solution, 
and  by  Cq  the  value  of  the  energy  integral  along  it. 

We  ask  ourselves  for  what  corrections  AXq,  AyQ,  AXq,  a£q  on  the 
initial  conditions  a  displacement  of  the  orbit  (3)  will  be  a  periodic 
orbit  with  about  the  same  period.  Let  us  denote  by  Tq  +  ATq  the  period 
of  this  varied  orbit,  and  by  Cq  +  ACq  its  constant  of  energy.  On 
expanding  the  periodicity  conditions 

x(To+iTo’xo+4xo’yo+A)'o,!io+Aio'yo+4^o)  ■  xo  +  too- 

y<VAVxo+Axo  ■  yo  +  Ayo- 

VTo+ATo-xo+Axo>yo+AVV  ■  xo  +  4xo’ 

yo  (To+ATo ,xo+Axo  ’  yo+iyo  ,xo+Axo  ,yo+Ayo)  ■  y0  +  Ay0* 


in  power  series  of  the  corrections  and  on  retaining  only  the  first  order 
terms,  we  come  to  the  linear  system 


x(Tq)-1 


y<V4x 

3^  i(T0)4x 


0 


0 


%  >(V 


0 


to0  +  3^  x<I0)Ay0  +  3^  x(VS  +  3^  x(T0)AyO  + 

+  y(T0)-1]Ay0  +  3i^  y(T0)iX0  +  3^  y(T0)Ay0  + 

+  3y^  X^VAy0  +  X^'V~1]  +  3y^  x^T0*AyO  + 

+  3^  *<VAy0  +  4  *<T0)4i°  +  [W0  ^V'1]  *»0  + 


XqAT  *  0, 
y0AT  -  0, 

XqAT  -  0, 

y0AT  -  0, 


(A) 


to  which  we  add  the  equation 
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(5) 


expressing  that  the  corrections  on  the  initial  conditions  result  in 
a  first  order  correction  ACq  on  the  constant  of  energy. 


Let  us  assume  that  the  rank  of  the  matrix 
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(6) 


is  equal  to  4.  Then  on  applying  Poincard's  method  of  continuity 
(Siegel  1956),  one  can  show  that  the  corrections  on  the  initial  conditions 
and  the  period  are  analytic  functions  of  ACq  in  the  neighborhood  of  the 
initial  energy  constant  Cq.  In  other  words,  there  exist  an  interval  0 
around  Cq  and  5  analytic  functions 
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(7) 


(8) 

on  that  interval  6  such  that  the  solutions  of  (1)  having  the  initial 
conditions  (7)  are  periodic  with  the  period  T  as  given  by  (8). 

Such  a  one-parameter  manifold  of  orbits  is  what  Wintner  (1931)  calls 
a  natural  family;  a  periodic  orbit  which  belongs  to  a  natural  family  is 
called  singular  by  Whittaker  (1916). 

Since  the  dynamical  system  is  conservative,  the  conditions  (7) 
imply  that  there  exist  two  sequences  x^(t),  yk(t)  of  functions  which 
are  periodic  with  period  T  such  that  the  series 


X0  -  *0  +  X  xo  k^C-C(P  * 

k>l  * 

Yo  ■  yo  +  ^  yo  k(c'co)k> 
k>l  ’ 

*0  -  4o  +  1.  io,k(c-co)k- 

k>l 

Yo  -  y0  +  .1,  *o,k<c-co)k> 

k>l 


T  ■  To  +  &  To,k<c-co> 


X(t)  -  x(t)  +  £  xk(t)(C-CQ)k, 
k>l 

(9) 

Y(t)  -  y (t)  +  l  yk(t)(C-C0)k 
k>l 


represent  the  natural  family  (7(C)  in  the  neighborhood  of  its  element 
0(CQ)  given  by  (3). 
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A  natural  family  of  periodic  orbits  defines  In  the  phase  space  a 
two  dimensional  torus  upon  which,  given  a  convenient  definition  of  the 
initial  point  on  each  orbit,  the  time  t  and  the  Painlevd  constant  C 
constitute  a  system  of  analytic  coordinates. 

Solving  the  problem  of  continuing  analytically  the  torus  0(C) 
from  the  initial  orbit  0(Cq)  means  determining  the  series  (8)  and  the 
time  dependent  coefficients  x^  and  y^  in  the  series  (9).  In  a  few 
simple  cases,  these  functions  can  be  determined.  But  in  general  the 
solution,  in  this  most  complete  sense,  is  not  feasible.  To  start  with, 
the  generating  solution  0(Cq)  is  usually  obtained  by  integrating 
numerically  the  equations  of  motion  (1),  so  that  the  very  first  coefficients 
x(t)  and  y(t)  are  obtained  only  in  the  form  of  tables  where  only  a  finite 
number  of  points  along  0( Cq)  are  entered.  And  even  when  the  initial  orbit 
is  expressed  in  a  somewhat  more  explicit  way  as  a  function  of  the  time,  it 
is  most  often  not  possible  to  find  in  the  same  explicit  way  even  the  time 
functions  x^(t)  and  y^(t)  in  the  expansions  (9). 

Since  the  analytic  continuation  of  a  natural  family  is  generally 
intractable,  it  is  quite  important  to  settle  upon  methods  of  numerical 
continuation.  The  difference  between  analytical  continuation  and  numerical 
continuation  can  be  sensed  most  distinctly  in  geometric  terms. 

We  think  of  the  equations  (9)  as  defining  the  continuous  deformation 
of  a  periodic  orbit  on  the  torus  0(C).  Starting  from  the  location 

{x(t;CQ)  ,y(t;CQ)  ,x(t;CQ)  ,y(t;CQ)  :  0  <_  t  £  TQ) 

in  the  phase  space,  the  orbit  will  stretch  and  twist  under  the  deformation 
so  as  to  occupy,  when  the  constant  of  energy  reaches  the  value  C,  the 


f 
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location 

{X(t,C),Y(t,C),X(t,C),Y(t,C):  0  <  t  <  T} 

defined  by  the  analytic  expansions  (8)  and  (9)  in  the  phase  space. 
Continuing  analytically  the  natural  family  0(C)  means  following  this 
continuous  deformation  of  the  periodic  orbit  0(Cq)  taken  as  a  whole. 

In  the  course  of  the  deformation,  each  point  of  the  initial  orbit 
0( Cq)  follows  a  certain  path  on  the  torus  0(C).  It  is  defined  in  the 
phase  space  by  the  numerical  series  (7).  Once  this  path  is  known,  the 
periodic  orbits  can  be  determined  unambiguously.  The  methods  of 
numerical  continuation  propose  to  follow  the  displacement  on  the  torus 
of  only  one  point  of  the  generating  orbit.  Thus  they  can  be  regarded  as 


procedures  by  which  the  initial  conditions  and  period  of  the  periodic 
orbit  0(CQ)  are  corrected  to  yield  the  initial  conditions  and  period 
for  the  neighboring  orbits  in  the  family. 


Figure  1.  The  two  steps  of  a  numerical  continuation. 
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As  for  the  methods  of  analytical  continuation,  in  general  it  is 
not  possible  to  compute  all  the  numerical  coefficients  in  the  series 
(7)  and  (8).  The  methods  of  numerical  continuation  aim  rather  at 
determining,  or  even  less  at  estimating  approximately,  only  their  first 
few  coefficients,  hereby  determining  only  approximately  the  deformation 
path  of  the  initial  conditions  on  the  torus  (7(C) .  For  this  reason  they 
provide  basically  two  schemes,  one  for  predicting  initial  conditions  and 
a  period  ahead  of  a  known  periodic  orbit,  and  the  second  for  correcting 
the  extrapolation.  In  regard  to  a  natural  family,  both  predictor  and 
corrector  should  have  specific  properties. 


Extrapolated  values  for  the  initial  conditions  (x^,y^,x^,y^) 
and  the  period  at  the  energy  level  Cq  +  AC  should  be  at  least 

such  that,  to  the  first  order, 


+  4C  •  8C  W  > 


yl  '  y0  +  4C  '  3C  W  ’ 
*1  *  *0  +  4C  •  3C  W  ’ 


yl  ‘  y0  +  iC  ’  9C  W ' 


T1  “  T0  +  AC  ’  dc  W  * 


It  means  that  the  predictor  should  request  as  a  minimum  a  tangent  to  the 
natural  family  (7(C)  at  the  initial  position  on  the  original  orbit  (7(Cq) 
in  the  direction  of  increasing  energy  constants,  as  well  as  the  tangent  at 
the  point  Cq  to  the  period  curve  T(C)  belonging  to  the  family. 
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Evidently  we  could  compromise  and  replace  the  tangents  by  the  secants 
going  through  two  previously  evaluated  sets  of  initial  conditions  on  the 
torus.  But  what  is  gained  from  simplifying  the  predictor  in  this  way 
often  results  in  poor  convergence  during  the  corrector  part  of  the  scheme. 
Moreover,  in  approaching  an  element  of  the  family  which  is  a  branching 
orbit  on  the  analytic  manifold  0(C)  or  an  extremum  on  the  inverse  C(T) 
of  the  period  curve,  extrapolation  along  the  secants  is  apt  to  derail 
without  notice  the  representative  initial  conditions  from  their  deformation 
course  on  the  family  0(C)  onto  another  natural  family  in  the  neighborhood. 
Besides  as  we  shall  show,  there  is  no  essential  difficulty  in  extracting 
the  tangential  elements  needed  for  an  efficient  predictor  from  a  numerical 
integration  of  the  variational  equations. 

The  phase  state  obtained  by  shifting  the  initial  conditions  along  a 
tangent  to  the  torus  0(C)  by  an  amount  equal  to  AC  does  not  generally 
lie  on  the  natural  family.  At  any  rate,  it  lies  on  the  integral  manifold 
^(Cq+AC)  defined  by  the  energy  integral  (2),  and  if  AC  is  small  enough, 
it  lies  close  to  the  periodic  orbit  which  is  the  intersection  of  the  torus 
0(C)  and  the  integral  manifold  ^(CQ+AC) .  Therefore,  the  corrector  stage 
of  a  numerical  continuation  should  aim  at  bringing  this  estimate  of  the 
initial  conditions  closer  to  the  periodic  orbit,  without  causing  it  to 
leave  the  integral  manifold  H  to  which  it  should  belong.  On  the  other 
hand,  when  only  first  order  corrections  are  retained,  several  iterations 
might  be  necessary  to  correct  the  predicted  values;  therefore  they  should 
be  made  to  converge  in  quadratic  convergence.  At  last,  a  first  order 
displacement  of  the  initial  conditions  tangentially  to  the  orbit  (r) 
that  they  define  is  not  generally  going  to  bring  the  phase  state  closer 
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to  the  periodic  orbit  0(Cg+AC)  since  it  amounts  only  to  an  advance 
of  the  time  on  the  orbit  r.  Therefore,  one  should  use  exclusively 
corrections  by  which  the  phase  state  is  displaced  normally  to  the  orbit 
T  to  which  it  belongs.  We  show  here  how  corrections  having  these  three 
properties  (of  being  isoenergetic ,  quadratically  convergent  and  normal) 
can  be  computed  by  integrating  numerically  the  homogeneous  Hill's  equation 
for  isoenergetic  normal  displacements  as  well  as  the  accompanying 
quadrature  for  isoenergetic  tangential  displacements. 

2.  INTRINSIC  DISPLACEMENTS  OF  AN  ORBIT 

Let 

x  -  x(t),  y  -  y(t)  (10) 

be  a  given  solution  of  (1)  defined  over  the  time  interval  I.  It  is 
assumed  to  have  the  following  properties: 

a)  It  is  not  an  equilibrium  position; 

b)  The  trace  of  the  orbit  over  its  curve  of  zero  velocity  is 
empty. 

Under  these  assumptions,  the  norm  of  the  velocity  along  the  orbit 

V  -  (x^y2)*5  (11) 

is  nowhere  zero  for  as  long  as  t  is  in  the  interval  I,  and  the 
inclination  <p  of  the  velocity  vector  on  the  x-axis  is  defined 
unambiguously  by  the  formulae 


x  ■  V  cos 


y  ■  V  sin  <p . 


(12) 
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A  displacement  (or  variation)  of  the  solution  (10)  is  a  pair  of 
numerical  functions 


u  =  u(t),  v  =  v  (t) 


that  is  a  solution  of  the  linear  differential  system 


0  -  2Av  +  (W  +2A  y)u  +  (W  +2A  y)v, 
xx  xJ  xy  yJ 

V  »  -2Au  +  (W  -2A  x)u  +  (W  -2A  x)v; 

xy  x  yy  y 


(13) 


it  is  understood  that  the  function  A,  its  partial  derivatives  as  well 
as  those  of  W  are  to  be  expressed  in  (13)  by  means  of  (10)  as  functions 
of  the  time  along  the  orbit. 

The  variational  system  (13)  admits  the  integral 


Y 


Wu+Wv-xu-yv. 
x  y  J 


(14) 


The  normal  and  tangential  components  (n,p)  of  the  displacement 
(u,v)  are  defined  respectively  by  the  identities 


n  *  -u  sin  <p  +  v  cos  4> , 
p  ■  u  cos  4)  +  v  sin 


(15) 


THEOREM.  A  numerioal  function  n  ■  n(t)  is  a  normal  displacement  of 
(10)  belonging  to  the  integration  constant  y  in  (14)  if  and  only  if  it 
satisfies  the  linear  differential  equation 


ti  +  On  ■  2y  (A+  $)  /V 


(16) 


where 
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0  -  V/V  +  2(A+$)2  +  2A2  -  W  -  W  -  2V(A  sin  <J>-A  cos  <f>) .  (17) 

xx  yy  x  y 

In  wTztofc  case,  a  numerical  function  p  «  p(t)  is  the  corresponding 
tangential  displacement  of  (10)  for  the  same  integration  constant  y 
if  and  only  if  it  is  given  by  the  quadrature 

~  (p/V)  -  2(A+$)n/V  -  Y/V2.  (18) 

In  order  to  prove  the  theorem,  we  begin  by  differentiating  (1) 
with  respect  to  the  time,  so  as  to  obtain  that 


*xv  -  y'u  +  xv  -  yii  *  2A(xu+yv-xu-yv) 

+  (xv-yu)  (W  +W  +2A  y-2A  x)  . 

J  xx  yy  ir  y 

We  now  propose  to  reduce  (19)  to  the  equation  (16).  To  this  effect  we 
invert  the  linear  transformation  (15)  to  obtain  that 


(19) 


u  =  p  cos  (j)  -  n  sin  (j>, 
v  *  p  sin  (J)  +  n  cos 

Then  we  use  (12)  and  (20)  to  compute  that 


(20) 


xv  -  yu  »  Vn;  (21) 

we  differentiate  (12)  with  respect  to  the  time  and  we  use  (20)  to  obtain 


xu  +  yv 


Vp  +  Vn$>; 


(22) 


also  we  differentiate  (20)  with  respect  to  the  time  and  calculate  that 


xu  +  yv  ■  Vp  -  Vn$ . 


(23) 
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Next  we  differentiate  (12)  twice  with  respect  to  the  time  so  as  to 
produce 

’xv  -  *yu  -  (V-vi2)n  -  (2V$+V<|>)p;  (24) 


also  we  differentiate  (20)  twice  with  respect  to  the  time,  which  yields 


u 


•  • 
v 


••  ••  2  ** 
p  cos  $  -  ii  sin  4>  +  u$  -  2v$  -  v<J>, 

••  ••  *  2  *  *  v 

p  sin  +  n  cos  $  +  +  2u<j>  +  u4> , 


and  thus  we  are  able  to  compute  that 


iv  -  yii  ■  V(n-n$2+2£<j>+p<j>)  . 

Finally  by  substituting  (21),  (22),  (23),  (24)  and  (25)  in  (19),  we 
obtain  that 

ri  +  [V/V-2$(2A+$)-W  -W  -2A  y+2A  x]n  « 

xx  yy  x  y 

2(A+j>)(Vp-Vp)/V. 

There  remains  now  to  eliminate  p  and  p  from  (26).  On  using  (1) , 
(15)  and  (22),  we  check  that 


(25) 


(26) 


Wxu  +  Wyv  -  Vp  +  Vn<(>  +  2VAn.  (27) 

Hence  (23)  and  (27)  helps  writing  the  integral  (14)  in  the  form 

Vp  -  Vp  -  -y  +  2V(A+4»)n.  (28) 

Eventually  we  substitute  (28)  in  (26) ,  and  this  yields  the  equation  (16) 
as  announced  in  the  theorem.  Of  course,  (28)  is  nothing  but  another  form 
of  the  quadrature  (18) . 
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That,  from  any  particular  solution  of  (16)  and  any  resulting 
quadrature  (18) ,  one  is  able  to  construct  a  solution  of  the  displacement 
equations  (13)  satisfying  the  first  integral  (14)  is  an  elementary  point 
whose  proof  we  leave  to  the  reader. 

Concerning  the  above  theorem,  the  following  comments  might  be 
appropriate: 

1)  Jacobi's  equations  (13)  which  determine  the  displacements 
(u,v)  constitute  a  Lagrangian  system.  In  that  respect,  the  theorem  can 
be  interpreted  as  providing  a  method  by  which  the  non-conservative 
integral  (14)  is  used  in  order  to  reduce  the  Lagrangian  system  of  Jacobi's 
equation  (originally  of  order  4)  to  a  linear  differential  system  of  order 
2,  to  be  followed  by  a  quadrature.  Straightforward  differentiations  and 
eliminations  were  sufficient  here  to  serve  this  purpose,  mainly  for  the 
reason  that  the  equations  to  be  reduced  as  well  as  the  reducing  integral 
are  linear. 

2)  The  function  0  defined  by  formula  (17)  is  well  known  in  various 
problems  of  Celestial  Mechanics. 

For  instance,  when,  in  our  expresssion,  we  take  A  as  a  constant 

function,  then  A  =  0  and  A  =  0,  so  that 

x  y 

0  *  V/V  +  2(A+i)2  +  2A2  -  Wxx  -  wyy. 

This  expression  can  be  found  in  Plummer  (1911) . 

Also  on  using  the  differential  equations  (1)  to  compute  that 

2  2  *2  2  *2 

W  +  W  «  V  +  V  (2A+$)  , 
x  y 


I 
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we  deduce  that 


Then,  since 


-2 

x 


vV  +  v2 


xx  +  yy 


..2 
+  x 


we  arrive  at  the  identity 


V 

V 


•2+x^> 

V 


Therefrom  we  obtain  by  substitution  in  (17)  another  form  for  0,  namely 


0  «  (xx+yyO/V2  +  3$2  +  4A$  +  4A2  -  W  -  W  -  2V(A  sin  <J>+A  cos  <fr)  .  (29) 

xx  yy  x  y 

When  A  =  0,  it  restitutes  the  expression  found  by  Poincard  (1899)  while, 
for  A  being  a  constant  function,  it  is  the  form  originally  given  by  Hill 
for  his  Lunar  Theory  and  by  Message  (1959)  for  the  planar  Restricted  Problem 
of  Three  Bodies. 

Again  on  observing  that 


W2  -  W2  +  2A(W  y-xW  )  -  V2  +  V2$2  +  2VA^, 

x  y  *  y 


we  conclude  that  0  can  also  be  defined  by  the  formula 


V20  -  VV  -  2V2  +  2(W2+W2)  +  4A(W  y-W  x)  +  4A2 

x  y  x  y 


-W  -  W  -  2 (A  y-A  x) . 
xx  yy  x'  y 


(30) 
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the  problems  they  are  concerned  with  (namely  the  mean  motion  of  the 
lunar  perigree  or  the  stability  of  a  periodic  orbit)  refer  exclusively 
to  the  case  when  0  is  a  periodic  function  of  the  time,  and  the  only 
quantities  they  are  interested  in  are  the  characteristic  exponents  for 
the  monodromy  matrix,  which  is  the  matrizant  of  precisely  the  homogeneous 
equation  (31)  at  the  end  of  the  period.  Poincar£  and  Birkhoff  make  an 
exception  in  that  respect,  as  both  adopted  the  Maupertuisian  form  of  the 
Principle  of  Least  Action  as  the  premiss  from  which  one  ought  to  derive 
Hill's  homogeneous  equation  (31).  Thereby  they  were  necessarily  confined 
to  isoenergetlc  displacements,  and  from  our  theorem,  it  results  quite 
clearly  that  the  Principle  of  Least  Action  can  lead  only  to  Hill's 
homogeneous  equation.  As  for  the  nonhomogeneous  equation  (16),  it  can 
only  be  derived  from  Hamilton's  general  Principle  of  Variation  wherein  the 
varied  motion  is  not  subject  to  the  restriction  that  on  it  C  remains 
constant . 

In  Poincare's  derivation  of  Hill's  homogeneous  equation,  several 
unsatisfactory  points  were  amended  by  Wintner  (1930).  By  this  expedient, 
Wintner  came  to  recognize  that  a  direct  approach  could  be  substituted  to 
Poincare's  oblique  treatment.  However,  because  he  still  adheres  to  an 
algebraic  identity  established  by  Poincard,  Wintner  (1931a)  makes  his 
construction  unnecessarily  clumsy.  We  like  to  acknowledge  that  our 
theorem  is  the  offshoot  of  an  effort  aiming  at  simplifying  Wintner 's 
argumentation. 

A  straightforward  deduction  of  the  homogeneous  equation  (31)  has 
been  proposed  by  Darwin  (1897),  using  the  arc  length  along  the  orbit  as 
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A  similar  evaluation  of  0  by  Poincar&  in  the  case  when  A  is  not 
uniformly  zero  was  proved  to  be  wrong  by  Wintner  (1930) .  On  correcting 
PoincarA  on  this  point,  Wintner  produced  the  formula  (30)  first  for  the 
constant  function  A,  and  then  generalized  it  to  any  function  A  which 
is  sufficiently  differentiable  (Wintner  1931a). 

3.  ISOENERGETIC  DISPLACEMENTS 

Since  it  is  linear,  equation  (16)  has  to  be  solved  in  two  steps. 

In  the  first  step,  one  has  to  find  the  matrizant  of  the  homogeneous 
equation 


n  +  0n  *  0 

obtained  from  (16)  by  omitting  the  right-hand  member.  In  point  of  fact 
this  amounts  also  to  putting  y  =  0  in  the  right-hand  member  of  (16). 
Therefore,  the  solutions  of  (31)  with  the  functions  determined  from  the 
resulting  quadrature  (18)  are  the  normal  and  tangential  components  of 
displacements  which  belong  to  the  particular  value  y  *  0  in  the  integral. 
For  this  reason,  they  are  called  isoenergetic  by  Wintner  (1931a). 
Conversely,  of  course,  the  normal  component  of  an  isoenergetic  displacement 
is  a  solution  of  the  homogeneous  equation  (31). 

This  exceptional  character  of  the  isoenergetic  normal  displacements 
has  not  been  recognized  quite  distinctly  by  Hill  in  his  Lunar  Theory.  As 
a  consequence  later  authors  lost  it  out  of  sight,  and  several  of  them 
resorted  to  unconvincing  arguments  in  order  to  justify  the  selection  of 
isoenergetic  displacements.  Actually,  in  relation  to  the  equation  (16), 


(31) 
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the  Independent  variable.  Darwin's  treatment  has  been  simplified  by 
Plummer  (1918).  Professor  Danby  has  informed  us  that  he  has  extended 
Plummer's  argumentation  to  encompass  the  nonhomogeneous  equation  (16). 

As  we  know,  the  general  solution  of  (31)  is  a  linear  combination 

n(t)  *  anI(t)  +  6nII(t)  (32) 

of  two  particular  solutions  satisfying  respectively  the  initial 
conditions 


nX(0)  -  1,  nX(0)  -  0, 

nn(0)  «  0,  nH(0)  -  1. 

The  real  numbers  a  and  8  are  the  two  arbitary  constants  of 
integration.  Besides,  since  the  Wronskian  of  any  pair  of  solutions  for 
(31)  is  a  constant,  we  have  that 


(33) 

(34) 


nI(t)nII(t)  -  nI(t)nII(t)  ■  1. 
Therefore,  the  matrizant  of  (31)  which  is  the  matrix 


(35) 


/  nX(t)  nXI(t)\ 
N(t)  -  [  ] 

W(t)  nU(t)/ 


is  unimodular,  so  that  its  inverse  is  simply  the  matrix 


(36) 


N"1(t) 


'hU(t)  -nII(t)' 


‘-h1 (t)  nX(t) 
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Consequent  ly  from  the  expression  of  the  general  solution  of  (16), 
namely 

-  N(t)  /  °\  +  2yN(t)  f  N_1(s)(  .  W 

\n0/  J0  \  [A(s)+<j>(s) ] /V(s)y 

we  conclude  that  the  function 


n(t)  -  nI(t)nQ  +  nn(t)h0  +  2Y[nI(t)a(t)+nII(t)b(t)] 


(37) 


where 


a(t)  -  -  f  nIX(s)  [A(s)+4>(s)  ]  /  V(s)ds 
J0 

ft  i 

b(t)  *=  In  (s)  [A(s)+<j>(s)]/V(s)ds 
0 

is,  among  all  normal  displacements  defined  by  the  initial  conditions 
n^  and  n^,  the  one  which  gives  the  value  y  to  the  integral  (14). 


(38) 


In  most  instances,  neither  0  nor  the  right-hand  member  of  (16) 
can  be  produced  in  literal  form,  hence  the  basic  isoenergetic  normal 
displacements  n*  and  n^*  can  be  obtained  only  by  numerical  integration. 
As  a  result,  one  should  look  for  a  better  way  of  obtaining  the  general 
normal  displacement  (32)  than  by  performing  the  quadratures  (38) .  Let 
us  consider  the  equation 


ii  +  On  =  2(A+$)/V 


(39) 


and  denote  by 


III 

n 


its  solution  determined  by  the  initial  conditions 


nm(0) 


Am(0) 


0. 
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Then  yn111  is  the  particular  solution  of  (31)  which  satisfies  the 
initial  conditions  n(0)  *  h(0)  *  0,  and,  in  view  of  (37),  the 
general  solution  of  (16)  can  also  be  given  the  form 

n  «  an  +  gn  +  yn  (40) 

where  a  and  8  are  two  arbitrary  constants  of  integration. 

4.  ISOENERGETIC  CORRECTOR 

Let  us  assume  that  x^,  y^,  Xq,  y^  are  approximately  the  initial 
conditions  for  a  periodic  orbit  belonging  to  the  energy  constant  C^; 
let  T  denote  an  approximate  estimate  of  its  period.  We  represent  by 
x(t),  y(t)  the  solution  of  the  equations  (1)  which  satisfies  the  initial 
conditions  xQ,  yQ,  xQ  and  yQ. 

The  problem  is  to  determine  the  first  order  corrections  u(t),  v(t) 
such  that 

X(t)  =  x (t )  +  u(t),  Y(t)  *  y(t)  +  v ( t ) 

be,  for  the  same  Painleve  constant  C^,  a  closer  approximation  to  the 
periodic  orbit.  The  corrections  to  the  initial  orbit  imply  that  the 
approximate  period  T  will  be  modified  at  the  first  order  by  an  amount 
AT. 

Let  denote  the  inclination  on  the  x-axis  of  the  velocity  at 

the  initial  time  t  =  0  along  the  initial  orbit.  By  orthogonal  projections 
on  the  initial  tangent  and  the  initial  normal,  we  obtain  the  periodicity 
conditions  in  the  following  form: 
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X(T+AT)cos  <|>q  +  Y(T+AT)sin  ^  ■  X(0)cos  +  Y(0)sin 
-X(T+AT)sin  <J)q  +  Y(T+AT)cos  <f>0  *  -X(0)sin  <f>0  +  Y(0)cos  <J>q, 
X(T+AT)cos  +  Y(T+AT)sin  ^  =  X(0)cos  <{)q  +  Y(0)sin  , 
-X(T+AT)sin  4>q  +  Y(T+AT)cos  <j>g  =  -X(0)sin  +  Y(0)cos  0^. 

In  these  expressions,  we  perform  the  substitutions 


X(T+AT)  «  x(T)  +  x (T) AT  +  u(T), 
Y(T+AT)  »  y (T)  +  y(T)AT  +  v(T), 
X(T+AT)  -  x(T)  +  x(T)AT  +  u(T), 
Y(T+AT)  -  y (T)  +  y(T)AT  +  v(T). 


On  omitting  terms  of  order  higher  than  one,  we  arrive  at  the  correction 
equations 


[u(T)-u(0)  ]cos  4>q  +  [v(T)-v(0)  ]sin  $0  +  VQAT  -  , 

[u(T)-u(0)]sin  <pQ  -  (v(T)-v(O)  ]cos  <j>0  ■=  A2> 

[u(T)-u(0)]cos  $0  +  [v(T)-v(0)  ]  sin  4>Q  +  VQAT  =  A3> 

[u(T)-u(0)]sin  -  [^(T)-^(O)  ]cos  -  V^aT  =  A*, 

where  we  put 

=  -fx(T)-x0]cos  <PQ  -  [y (T)-yQ] sin  $Q, 

A2  -  -[x(T)-Xq] sin  <J>0  +  [y(T)-yQ]cos  $Q, 

A3  =  -[x(T)-xq]cos  -  [y(T)-yQ]sin  <J>0 , 

A^  «  “[x(T)-xq] sin  +  [y (T)-yQ] cos  $Q. 


(41) 


(42) 
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The  corrections  will  be  decomposed  into  their  normal  and  tangential 
components  along  the  initial  orbit  so  that,  at  the  first  order,  we  have 
that 

u(T)  -  u(0)  =  ip (T)-p(O) ]cos  4 iQ  -  [n(T)-n(0)]sin  4>Q, 

v (T)  -  v(0)  =  [p(T)-p(0)  ] sin  4jq  +  [n (T)-n(O)  ]cos  4>0 , 

u(T)  -  u(0)  =  [p(T)~p(0)  ]cos  4)q  -  [n(T)-ft(0)]sin  4>Q 

(43) 

-  4>q{  [ p (T)  — p  (0)  ] sin  4>0  +  [n(T)-n(0)]cos  4>Q } , 

v(T)  -  v(0)  =  [p(T)-p(0)]sin  4^  +  £n (T) -n (0)  ]cos  <f>Q 

+  4>0{  [p(T)-p(0)  ]cos  4>q  -  [n(T)-n(0)]sin  4>Q} . 

We  substitute  the  expressions  (43)  in  the  correction  equations  (41) 
and,  after  a  few  manipulations,  we  arrive  at  the  system 

n(T)  -  n(0)  =  -&2, 
n(T)  -  h(0)  =  -A^  - 

P (T)  -  p(0)  +  VqAT  =  Ax,  (44) 

P (T)  -  p(0)  +  VqAT  =  A3  -  A2iQ. 

We  should  now  observe  that  these  four  correction  equations  are  not 
independent.  Indeed  from  the  variational  integral  in  the  form  (28),  we 
deduce  immediately  that 

V0[p(T)-p(0)]  =  V0[p(T)-p(0)]  +  2V0(A0+4>0)[n(T)-n(0)], 

which  proves  that,  at  the  first  order,  the  fourth  of  the  equations  (44) 
ought  to  be  a  linear  combination  of  the  first  and  the  third  one.  This 
implies  that,  at  the  first  order,  we  ought  to  have 
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Vi  -  V0(2A0^0)A2  ‘  V3  =  °*  (45) 

Thus  by  computing  this  quantity  at  the  period  T  along  the  initial 
orbit,  we  shall  gain  an  estimate  of  the  second  order  contributions  to 
the  corrections,  which  contributions  we  have  decided  to  omit  from  the 
corrector  scheme.  And  it  is  one  of  the  aims  of  the  iteration  on  the 
corrections  to  bring  these  contributions  down  as  close  to  zero  as 
possible. 

The  variational  integral  (28)  informs  us  also  that  for  the 
corrections  on  the  initial  orbit  to  be  isoenergetic  to  the  first  order, 
it  is  necessary  and  sufficient  that  the  displacements  n  and  p  be  in 
turn  isoenergetic.  As  we  settled  down  for  an  isoenergetic  corrector, 
from  now  on  in  this  paragraph,  we  shall  restrict  ourselves  to 
isoenergetic  displacements. 

Thus  the  normal  correction  will  present  itself  as  the  linear 
combination 

I,  II 

n  =  c»n  +  Bn  .  (46) 

Consequently,  since  the  quadrature  (18)  is  linear  in  p  and  n,  the 
tangential  correction  will  be  of  the  form 

p  =  ap1  +  Bp11  (47) 

where  p*  (resp.  p1*)  results  from  n*  (resp.  n**)  in  virtue  of  the 
quadrature 

£<v>  -  2<A+»>?  • 


(48) 
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Before  we  fix  the  initial  conditions  p*(0)  and  p^^(0)  requested 
by  the  quadratures,  we  should  realize  that,  since  the  dynamical  system 
is  conservative,  the  selection  of  an  initial  value  p(0)  for  the 
tangential  displacement  amounts  to  nothing  else  than  a  translation  of 
the  time  origin,  which  is  an  irrelevant  operation  when  the  problem  is 
to  determine  the  initial  conditions  and  the  period  of  a  periodic  orbit. 
Therefore,  without  any  loss  of  consistency,  we  can  impose  that 

P^O)  =  pn(0)  =  0.  (49) 

In  geometrical  terms,  the  only  displacement  that  our  corrector  provides 
for  the  position  is  a  shift  along  the  normal  to  the  original  orbit  at 
the  initial  position. 

On  using  the  initial  conditions  (33),  (34),  (49)  and  on  substituting 
(46)  and  (47)  in  the  correction  equations  (44),  we  readily  obtain  that 
the  corrective  displacement  a  along  the  normal  and  the  modification  B 
for  the  initial  normal  velocity  are  solutions  of  the  linear  system 

a  [nI(T)-l]  +  Bnn(T)  =  -A^ 

anI(T)  +  8[nn(T)-l]  -  -A4  -  A^.  (50) 

In  view  of  (35) ,  its  determinant  is  found  to  be  equal  to 

2  -  [nI(T)+nII(T)J . 

Thus  the  cases  when  the  trace 


Tr(T)  =  n1  (T)  +  nH(T) 
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of  the  matrizant  N(T)  for  Hill's  homogeneous  equation  is  equal  to  2, 
turn  out  to  be  singularities  for  our  corrector  scheme.  As  we  rule  them 
out  of  order  here,  we  may  assume  that  the  system  (50)  has  a  unique 
solution.  Then  much  in  the  same  way  as  we  did  for  the  system  (50),  we 
extract  from  the  third  of  the  equations  (44)  the  first  order  correction 
to  the  period: 

AT  -  [A1-ctpI(T)-8pII(T)]/V0. 

Once  the  factors  a  and  8  have  been  determined,  we  can  compute 
the  corrections  to  the  initial  conditions  in  Cartesian  coordinates  through 
the  following  sequence  of  formulae 

Axq  =  u(0)  =  -a  sin  $Q, 

AyQ  =  v(0)  =  a  cos  <pQ, 

Axq  =  u(0)  =  a[2A(x0,y0)+^0]cos  <j>Q-  8  sin  <frQ, 

AyQ  ■  v(0)  *  a[2A(x0,y0)+$0]sin  +  8  cos  <t>Q. 

On  inserting  the  new  initial  positions 

xo  +  Axo’  yo  +  Ay0 

into  Painlev^'s  integral  (2)  where  the  constant  C  is  given  the  value 
Cq,  we  compute  the  new  norm 

v0  ■  l2w(xo+ix0'y0+ay0)"c01'i 

of  the  velocity  vector.  But  the  quantities 

*0  +  S’  y0  +  iy0 


(51) 


) 
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give  the  new  direction  of  this  vector.  In  consequence,  we  compute  the 
norm  of  this  direction  vector 


V0  -  («0+Ax0)2  +  (yo+ayo)2]*5 


so  that  for  the  new  initial  velocities  we  obtain  the  numbers 


v.  idi)  v,  fa^.o 


0  V, 


0  V, 


5.  TANGENTIAL  PREDICTOR 


Given  a  periodic  orbit 


^WVVW-  y(t;to,WVyo!Co) 


with  period  T^,  which  is  determined  at  time  t^  by  the  initial 
conditions  suc^  t^iat  Painleve's  constant  takes  the 

value  CQ,  let  us  assume  that,  for  this  orbit,  the  matrix  (6)  is 
of  rank  4.  Thus  there  exists  in  the  neighborhood  of  Cq  a  natural 
family  0(C)  of  periodic  orbits  which  contains  the  given  periodic 
orbit.  Moreover,  since  the  continuation  0(C)  is  analytic,  the  family 
is  represented  by  expansions  of  the  form 


X(t;CQ+AC)  =  x(t  ;Cq)  +  AC  •  x(t;CQ)  +' 


Y(t;C0+AC)  =  y(t;CQ)  +  AC~y(t;C0)  +■ 


(53) 


and  the  period  T  along  the  natural  family  likewise  takes  the  form  of 


a  series  in  AC: 
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T(C0+AC)  -  T0  +  AO^T0  +•••  .  (54) 

The  first  order  of  the  series  (53)  have  to  check  to  the  first  order  the 
variational  integral 

Wx(x(t:C0),y(t:C0>)  1c  x(t;C0)  +  Wy<*<t:C0),y(t;C0))  ^y(tjc0) 

•  x(t’C0)  X  i(t:C0)  •  #<t!C0)  sc  '  v 

Basically  we  do  concern  ourselves  here  with  the  problem  of  actually 
computing  as  function  of  the  time  the  displacement  (3x(t ;Cq) /9C, 9y (t ;Cq) /9C) 
which  causes  the  varied  orbit  to  be  periodic  with  the  modified  period  (54). 

For  the  sake  of  convenience,  we  rather  put 

y  -  >54C,  »  -  2  f|,  V  ’  2  !c  >  1T  '  2  Sc’  ■ 

In  these  notations,  the  problem  can  be  reformulated  as  that  of  finding  a 
displacement  (u,v)  which  satisfies  the  variational  integral 

Wu+Wv-xu-yv=l  (55) 

x  y  J 

and  Is  such  that,  for  y  sufficiently  small,  the  orbit 

X(t)  =  x(t)  +  yu(t)  +•••  Y(t)  =  y(t)  +  yv(t)  +••• 
be  periodic  with  period 


T  -  Tq  +  yil. 

Exactly  as  we  did  when  we  set  up  the  corrector  scheme,  we  project  the 
first  order  conditions  of  periodicity  orthogonally  onto  the  tangent  and 
the  normal  at  the  initial  point  of  the  given  orbit  O(Cq)  •  Introducing 
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therein  the  tangential  and  normal  components  of  the  displacement,  we 
eventually  arrive  at  the  following  conditions  of  periodicity 

n(T)  -  n(0)  =  0, 
n(T)  -  n(0)  «  0, 
p(T)  -  p(0)  +  VqAT  =  0, 

P (T)  -  p(0)  +  VqAT  -  0. 

But,  in  view  of  (55),  the  normal  and  tangential  displacements  satisfy 
the  first  integral 

Vp  -  Vp  -  2V(A+$)n  -  1. 

It  implies  in  particular  that 

V0[p(T)-p(0)]  =  Vo[p(T)-p(0)J  +  2V0(A0+<^0)[n(T)-n(0)], 

so  that  the  third  condition  of  periodicity  turns  out  to  be  a  linear 
combination  of  the  first  and  the  second  one. 

Now  from  (57),  we  also  conclude  that  the  normal  displacement  is  a 
solution  of  Hill’s  nonhomogeneous  equation 

ri  +  0n  *  2 (A+4>) /  V.  (58) 

Then,  as  we  have  shown  above,  it  is  a  linear  combination  of  the  type 

I  .  II  .  Ill 
n  *  an  +  gn  +  n  , 

where  n*  and  n**  are  the  basic  solution  of  Hill’s  homogeneous 
equation,  and  n**^  is  the  particular  solution  of  (58)  for  the  initial 


(56) 


(57) 
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conditions  n^*(0)  *  n*XX(0)  ■  0.  Correspondingly  the  tangential 
displacement  is  the  linear  combination 


1  ,  .  n  .  Ill 
p  *  ap  +  6p  +  p 


where  p11*  is  given  by  the  quadrature 


_d_ 

dt 


,  A+4  III  1 
2  -rp  n  -  — 

V 


for  the  usual  initial  conditon 


0. 


As  a  result  of  all  preceding  remarks,  the  conditions  of  periodicity 
(56)  becomes 

a[nX(T)-l]  +  Bnn(T)  =  -nin(T), 
anI(T)  +  6[nII(T)-l]  =  -nIIX( T) , 
apI(T)  +  8pn(T)  +  VqAT  =  -pHI(T). 


The  determinant  of  the  first  two  equations  is  equal  to  2  -  Tr(T). 

For  the  sake  of  brevity  we  shall  not  discuss  here  the  singular  case  of  a 
periodic  orbit  for  which  Tr(T)  *  2,  i.e.  of  an  orbit  whose  characteristic 
exponent  is  equal  to  zero.  Let  it  be  mentioned  that  it  occurs  only  in 
exceptional  circumstances,  like  a  relative  extremum  for  the  period  function 
T(C),  a  bifurcation  on  the  natural  family  0(C)  or  even  an  essential 
singularity  beyond  which  no  analytical  continuation  of  the  manifold  is 
permissible,  which  is  what  Wintner  (1931a)  calls  a  natural  termination  of 
the  family.  In  principle  these  situations  have  been  analyzed  by  Poincard 
(1899).  Recently  they  have  received  more  careful  attention  from  Hdnon 
(1965)  in  relation  to  several  natural  families  of  synodically  symmetric 


(59) 


orbits  in  the  Restricted  Problem  of  Three  Bodies. 
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Once  the  coefficients  a  and  8  have  been  obtained,  we  can  compute 
the  derivations  at  C  =  Cq  of  the  initial  conditions  which  determine 
the  periodic  orbit  0(Cq) »  namely 


Tc  W  *  2  u(0) 


“  T  0  sin 


W  ’  2  V<0)  ‘  2  “  cos  V 


dC  W  *  2  “<0)  * 


1  cos  *0  1  •  1 
2  ~y~  +  y(2A0+$0)a  C0B  *0  "  2  6  sin  V 


to  W  ■  2  *(0)  ’ 


1  s*n  1  1 
2  — y - +  2'(2Ao‘hVa  sin  ^0  +  ~2  3  COS  * 


0‘ 


In  geometric  terms,  we  hold  here  a  tangent  to  the  manifold  0(C)  at 
the  point  C^.  Thus  in  order  to  obtain  a  first  order  approximation  to 
the  next  orbit  in  the  family  which  lies  on  the  manifold  ^(Cq+  O 
determined  by  the  integral  (2),  we  suggest  moving  the  initial  point  to 
the  point  whose  coordinates  are 


xo  =  xo  +  4C-dC  xo(co)> 


y0  =  y0  +  40  '  dC  W  ‘ 


Then  we  compute  the  norm  of  the  velocity  vector  from  Painleve's  integral 

V*  -  [2W<x*,y*)-C0-AC]%. 

At  last  we  estimate  the  direction  of  this  vector  by  the  two  numbers 
X0  =  *0  +  AC  *  dC  *0^0^  *0  "  *0  +  AC  *  dC  * 
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so  that,  if  we  put 


V, 


(x2  +  y2)^ 
v  0  y0;  ’ 


the  new  initial  velocity  will  have  for  components 


X0  yo 

x*  =  V*  -zr-  v*  =  V*  — 

X0  V0  vn  y0  Oy 
u  0 


At  this  stage,  then  we  resort  to  the  corrector  scheme  to  move  the 
newly  obtained  initial  conditions  (point  B  in  Fig.  1)  by  successive 
approximations  along  the  manifold  iJ( Cg+AC)  right  onto  those  of  the 
periodic  orbit  0(Cq+AC)  (point  A'  in  Fig.  1). 


Figure  2.  Diagram  of  numerical  continuation  for  a  natural  family  of 
periodic  orbits. 


w  H,'  ?  V 
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6.  CHARACTERISTIC  EXPONENTS 

As  a  by-product  of  our  numerical  computation,  we  collect  immediately 
the  nontrivial  characteristic  exponents  of  each  periodic  orbit  that  we 
compute. 

Let  to  be  one  of  these  exponents.  There  exists  a  displacement  of 
the  orbit  which  is  of  the  form 

u(t)  -  eutU(t),  v(t)  =  eU)tV(t)  (61) 

where  U  and  V  are  periodic  functions  with  period  T.  We  prove  that 
the  displacement  (61)  is  isoenergetic .  Indeed,  along  the  varied  orbit, 
the  first  order  variation  of  the  energy  integral  is  of  the  form 

Y  *  P(t)ewt  (62) 

where  P  is  a  periodic  function  with  period  T.  Since  this  variation 
is  an  integral  of  the  variational  equations,  (13)  implies  that  P  is 
identically  equal  to  zero,  which  means  that  the  constant  y  is  equal  to 
zero,  and  this  establishes  that  (61)  is  an  isoenergetic  displacement. 

Now  the  normal  component  of  (61)  is  of  the  form 

n(t)  -  ewtN(t)  (63) 

where  N  is  a  periodic  function  with  period  T.  But  (63)  expresses 
precisely  that  is  a  characteristic  exponent  of  Hill's  homogeneous 
equation. 

It  follows  that,  in  order  to  compute  the  nontrivial  characteristic 


exponents  of  a  periodic  orbit,  it  is  necessary  and  sufficient  to 
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calculate  the  characteristic  exponents  of  Hill’s  equation  associated 
with  it. 

But  as  we  know,  these  quantities  are  such  that 

ujT  ,  -u)T 

s^  =  e  and  =  e 

are  the  eigenvalues  of  the  matrizant  N(T)  given  by  (36)  at  the  end  of 
the  period  of  the  function  0,  which  is  precisely  the  period  of  the  orbit 
Therefore  s^  and  are  the  roots  of  the  characteristic  equation 

s2  -  [nI(T)+nII(T)]s  +  [nI(T)nII(T)-nII(T)nI(T)  ]  =  0. 

In  view  of  the  Wronskian  integral  (35) ,  it  can  even  be  written  simply  as 

s2  -  Tr (T)s  +1=0. 

Consequently,  the  trace  Tr(T)  can  be  used  as  an  index  of  the 
stability  of  the  orbit. 

a)  If  j Tr (T) |  >  2,  the  characteristic  exponents  of  the  orbit 
are  of  the  unstable  type. 

b)  If  | Tr (T) |  <  2,  they  are  of  the  stable  type. 

c)  Tr(T)  =  +2  or  -2  represents  what  Wintner  terms  the  two 
cases  of  indifferent  stability. 


7.  SYMMETRIC  ORBITS 

Let  us  assume  that  the  functions  W  and  A  which  characterize  the 
dynamical  system  have  the  following  symmetry  property: 


(64) 


W(x,-y)  =  W(x,y), 


A(x,-y)  -  A(x,y) 
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/ 


identically  in  the  configuration  plane. 

Then  the  substitutions 

t  -*•  -t,  y  -*•  -y 

leave  the  equations  of  motion  unchanged.  This  invariance  in  turn  implies 
that,  if  the  pair  x(t),y(t)  is  a  solution  of  the  equations  of  motion 
(1),  then  the  pair  £(t),n(t)  such  that 

£(t)  =  x(-t),  n(t)  =  -y(-t) 


is  also  a  solution  ("Principle  of  Symmetry").  In  particular,  a  solution 
defined  by  the  initial  conditions  (xo*yo'^0*^0^  such  that 


y 


0 


(65) 


is  symmetric  with  respect  to  the  x-axis.  Indeed,  in  view  of  Cauchy's 
theorem  of  uniqueness,  the  Principle  of  Symmetry  results  here  in  the 
identities 


x(-t)  =  x(t),  y(-t)  =  -y (t) .  (66) 

It  is  easily  checked  that  the  conditions  (66)  which  define  an  orbit 
symmetric  with  respect  to  the  x-axis  imply  the  following  identities 

V(-t)  -  V(t),  cos  4) (— t )  *  -cos  4> (t )  , 

$>(-t)  -  4>(t),  sin  4>(-t)  *  sin  4>(t). 

As  a  result,  along  such  a  symmetric  orbit,  the  coefficient  0  in  Hill's 
equation  is  an  even  function  of  the  time. 
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Now  we  propose  to  show  that,  along  the  symmetric  orbit  defined 
by  the  initial  conditions  (65),  the  isoenergetic  normal  displacement 
n*  (reap,  n*1)  is  an  even  function  ( resp .  an  odd  function)  of  the 
time. 

Indeed  the  resolvent  N(t;0)  of  Hill’s  equation  is  characterized 
as  the  solution,  defined  by  the  initial  condition 

N(0;0)  =  I2, 

to  the  differential  matrix  equation 

N(t;0)  =  H(t)N(t ;0)  (67) 

where 

H(t)  =/  (68) 

'0(t)  0/ 


Since  0(-t)  *  0(t),  the  matrix  function  H(t)  is  an  even  function  of 
the  time.  Therefore,  the  substitution  t  ->  -t  in  (67)  provides  the 
identity 


N(-t;0)  =  -H(t)N(-t;0).  (69) 

If  we  put 


and  observe  that 


S 


2 


SH (t)S  =  -H(t)  , 
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we  see  that  (69)  can  also  be  written  in  the  form 

~  [ SN (— t ; 0) ]  =  H(t) [ SN(-t ;0) ] . 

Hence,  on  using  the  fact  that  N(t;0)  is  the  resolvent  of  (67),  we 
come  at  last  to  the  identity 

SN(-t;0)  =  N(t ;0)S.  (70) 

But  the  matrix  identity  (70)  is  equivalent  to  the  scalar  identities 

nX(-t)  -  nr(t),  nU(-t)  -  (71) 

This  concludes  the  proof  of  the  proposition. 

From  (71)  together  with  (18),  it  results  that 

pVt)  -  -p^t),  pIX(-t)  *  pH(t).  (72) 

Now  that  we  have  reviewed  the  symmetry  properties  of  the  displacements 
of  a  symmetric  orbit,  we  impose  the  condition  that  the  symmetric  orbit  be 
also  periodic.  Let  T  be  its  period. 
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Because  the  matrix  function  H(t)  is  periodic  with  the  same 
period  T,  on  substituting  t  +  T  for  t  in  (67),  we  obtain  that 

N(t+T;0)  =  H(t)N(t+T;0) 


and  deduce  therefrom  that 

N(T+t;0)  =  N(t;0)N(T;0) 

at  any  time  t.  In  particular,  for  t  =  -T/2, 

N(T;0)  *  N_1(-T/2;0)N(T/2;0). 

(73) 

(74) 

(75) 


This  matrix  formula  gives  rise  to  the  scalar  expressions 


nI(T)  =  hII(T)  =  nI(T/2)hII(T/2)  +  nII(T/2)nI(T/2) , 


nI(T)  =  2nI(T/2)nI(T/2), 


nIX(T)  »  2nII(T/2)nII(T/2) 


Note  that  the  Wronskian  integral  taken  at  half  the  period, 


nI(T/2)nII(T/2)  -  nI(T/2)nII(T/2)  =  1, 


) 
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gives  for  (73)  an  expression 

nX(T)  =  nIX(T)  =  1  -  2nI(T/2) 11 (T/2) 
that  will  be  useful  later  on. 

Formula  (73)  is  especially  interesting  for  the  computation  of  the 
characteristic  exponent  of  a  symmetric  periodic  orbit.  Indeed,  as  we 
have  seen  before, 

2  cosh  <jj  T  =  Tr (T)  =  nI(T)  +  nn(T), 
so  that,  in  view  of  (73), 

cosh  (i)  T  =  n1  (T/2) h11  (T/2)  +  n1  (T/2)nU (T/2)  . 

Thus  integrating  the  orbital  equations  and  Hill’s  homogeneous 
equation  for  only  half  a  period  is  sufficient  to  compute  the 
characteristic  exponent  u  of  a  symmetric  periodic  orbit. 

This  proposition  which  was  known  to  Darwin  (1897)  has  been 
rediscovered  by  Moulton  (1914)  and  applied  extensively  by  Lemaitre  in 
his  explorations  of  the  symmetric  periodic  solutions  to  Stdrmer's 
problem.  We  plan  to  show  that,  like  for  the  computation  of  the 
characteristic  exponent,  the  numerical  continuation  of  a  natural  family 
consisting  of  symmetric  orbits  requires  the  orbit  and  its  intrinsic 
variations  to  be  calculated  over  only  half  the  period. 

In  order  to  do  so,  we  need  to  back  up  and  consider  the  problem  of 
converging  iteratively  toward  the  initial  conditions  determining  a 


(76) 


(77) 
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non-necessarily  symmetric  periodic  orbit.  Here  we  shall  start  from 
the  periodicity  conditions 

X(T/2+AT/2)  «  X(-T/2-AT/2) , 

Y(T/2+AT/2)  =  Y(-T/2-AT/2) , 

X(T/2+AT/2)  =  X(-T/2-AT/2), 

Y(T/2+AT/2)  =  Y (-T/2-AT/2) . 

To  the  first  order  in  the  variations,  they  generate  the  correction 
equations : 

u(T/2)  -  u(-T/2)  +  [x(T/2)+x(-T/2)  ]AT/2  »  -[x(T/2)-x(-T/2) ] , 

v(T/2)  -  v(-T/2)  +  [y(T/2)+y(-T/2) ] AT/2  -  -[y(T/2)-y (-T/2) ] , 

u(T/2)  -  u(-T/2)  +  [X(T/2)-tfc(”T/2)  ]  AT/2  -  -[x(T/2)-x(-T/2) ] , 

v(T/2)  -  v(-T/2)  +  [y (T/2)+y (-T/2)  ]  AT/2  =  -[y(T/2)-y (-T/2) ] . 


By  orthogonal  projection  onto  the  tangent  and  the  normal  to  the 
orbit  at  the  time  T/2,  we  obtain  from  them  that 

p(T/2)  -  p(-T/2)  +  V(T/2) AT  =  A^ 
n(T/2)  -  n(-T/2)  =  -A2> 

(79) 

p(T/2)  -  p(-T/2)  +  V(T/2) AT  =  A3  -  A2*(T/2), 
h(T/2)  -  h(-T/2)  =  -A4  -  A1i(T/2), 


where  we  have  put 
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*  -[x(T/2)-x(-T/2)]cos  4»(T/2)  -  [y (T/2)-y (-T/2) ]sin  <j>(T/2), 

A9  =  -[x(T/2)-x(-T/2)]sin  $(1/2)  +  [y(T/2)-y(-T/2)]cos  $(T/2), 

(80) 

A3  =  -[x(T/2)-x(-T/2)]cos  4>(T/2)  -  [y (T/2)-y (-T/2) ]sin  $(T/2), 

A4  =  -[x(T/2)-x(-T/2) ] sin  4>(T/2)  +  [y(T/2)-y(-T/2)]cos  4>(T/2)  . 

If  we  use  the  correction  equations  (79)  to  set  up  an  isoenergetic 
corrector,  we  have  to  assume  that 

I  .  II  .  I  „  II 

a  =  an  +  6n  and  p  =  ctp  +  3p  , 

in  which  case  the  formulae  (79)  become 

a[nI(T/2)-nI(-T/2)]  +  B[nn(T/2)-nn(~T/2)  ]  «  -Aj, 

ot[nI(T/2)-hI(-T/2)]  +  B[nIX (T/2)-nn(-T/2) ]  -  -A4  -  A^d/2), 

II  II  II  (81) 

a[pi(T/2)-pi(-T/2)]  +  B[p  (T/2)-pXi(-T/2) ]  +  V(T/2)AT  -  A^ 

a[pI(T/2)-pI(-T/2)]  +  B[pII(T/2)-pII(-T/2) ]  +  V(T/2)AT  -  A3  -  A2$(T/2). 

Let  us  now  assume  that  the  generating  orbit  x(t),y(t)  is 
symmetric  with  respect  to  the  axis  Ox,  its  initial  conditions 
satisfying  the  relations  (65).  Then,  in  view  of  the  symmetry  identities 
(66),  we  obtain  that,  to  the  first  order  in  the  displacements, 

cos  $(T/2)  ■  0 

which  implies  that 


A^  -  -2y(T/2)sin  <j>(T/2),  A2  ■  0, 
A4  »  -2x(T/2)sin  $(T/2),  A3  -  0. 
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Taking  also  into  account  the  symmetry  identities  (71)  and  (72),  we 
are  thus  able  to  reduce  the  correction  equations  (81)  to  the  simpler 
system 

BnH(T/2)  =  0, 

anI(T/2)  =  (x(T/2)+y(T/2){(T/2)]sin  «*> (T/2)  , 
apI(T/2)  +  j  V (T/2) AT  =  -y(T/2)sin  <fr(T/2) , 

6pH(T/2)  =  0. 

This  system  needs  to  be  discussed.  If  n**(T/2)  =  0,  then  B 
is  left  undetermined  by  the  first  equation.  However,  in  view  of  (75) 
and  (76),  we  have  in  this  case  that  n^^(T)  =  0  and  n*(T)  =  1,  so 
that  we  find  ourselves  here  in  the  critical  case  where  Tr(T)  *  2. 
Moreover,  from  the  variational  integral  (28)  along  an  isoenergetic 
displacement  (y  =  0) ,  we  also  have  that  pII(T/2)  =  0,  since,  to  the 
first  order,  V(T/2)  =  0.  Thus  the  fourth  equation  in  (82)  is 
compatible.  Conversely,  if  p**(T/2)  =  0,  then  we  can  show  in  the  same 
way  that  n  (T/2)  =  0.  Hence,  if  we  exclude  the  critical  case  when 
Tr(T)  =2,  we  must  have  that  6=0. 

The  second  of  the  equations  (82)  has  a  unique  solution  if  and  only 
if  nI(T/2)  i  0.  But,  in  view  of  (74)  and  (76)  nI(T/2)  =  0  implies 
that  nJ(T)  =  0  and  consequently  that  Tr(T)  =  2.  Therefore,  excluding 
once  more  the  critical  case,  we  can  extract  from  the  second  of  the 
equations  (82)  the  correction  displacement  a  and  thereafter  from  the 
third  equation  the  correction  on  the  period.  It  is  then  obvious  that 


(82) 


the  corrections  on  the  initial  conditions  are 
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Axq  =  -a  sin  4>Q, 
Ay0  -  0, 


AyQ  =  a  [2A(Xq,0)  +  4>q] sin  <*>Q. 

Actually  the  last  relation  is  there  only  to  determine  the  sign  of  the 
correction  to  be  brought  to  the  initial  velocity  y^.  Indeed,  the  new 
velocity  orthogonal  to  the  axis  of  symmetry  should  be  computed  in  norm 
rather  from  the  energy  integral  to  make  sure  that  the  corrector  absorbs 
even  the  second  order  variation  of  the  energy  constant. 

Concerning  the  above  isoenergetic  corrector  for  symmetric  periodic 
orbits,  the  following  remarks  are  of  relevance. 

a)  The  fact  that  6  is  equal  to  zero  whenever  Tr(T)  t  2, 
implies  that  the  improved  orbit  generated  by  the  corrector 
from  a  symmetric  orbit  is  also  a  symmetric  orbit. 

b)  Here  no  use  is  made  of  the  fact  that,  at  half  the  period, 
the  orbit  ought  to  cross  orthogonally  the  axis  of  symmetry.  We 
know  from  experience  that,  unless  the  orbit  is  of  a  simple  shape, 
such  a  condition  makes  an  uneasy  criterion  of  periodicity.  For 
instance,  if  the  second  crossing  occurs  in  the  neighborhood  of 

a  singularity  or  close  to  a  curve  of  zero  velocity,  the  periodic 
orbit  will  show  intrincate  loops,  which  cross  several  times  the 
axis  of  symmetry  at  points  which  are  close  to  each  other  and 
with  inclinations  that  cannot  be  discriminated  distinctly.  In 
these  circumstances,  it  is  altogether  too  tempting  to  make 
subjective  decisions,  thereby  running  the  risk  of  jumping  into 


(83) 
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a  natural  family  different  from  the  one  to  be  continued.  In 
our  corrector,  we  converge  to  the  symmetric  periodic  orbit 
nimply  by  shifting  conveniently  the  initial  position  along 
the  axis  of  symmetry,  and  by  adjusting  the  period  to  this 
translation.  The  fact  that,  ultimately  when  the  procedure  has 
converged,  the  periodic  orbit  crosses  orthogonally  the  axis  of 
symmetry  is  not  a  boundary  condition  that  we  impose  to  decide 
the  time  at  which  we  stop  the  integration,  but  a  consequence 
of  an  iteratively  convergent  adjustment  of  the  initial 
conditions. 

When  we  use  the  correction  equations  (79)  for  a  tangential  predictor, 
we  have  to  assume  that 


n  = 


I  _ 
an  + 


„  11  ^  111 
Bn  +  n  , 


I  , 
ap  + 


0  II  .  Ill 
Bp  +  P 


and  that 


A1  *  A2  ‘  fi3  =  \ 


Thus  we  arrive  at  the  formulae 

a[nI(T/2)-nI(-T/2)J  +  B[nII(T/2)-nII(-T/2)  ]  =  -[nm(T/2)-nm(-T/2)  ] , 
a[hI(T/2)-hI(-T/2)]  +  B[nII(T/2)-hII(-T/2)]  =  -[h111 (T/2)-n111  (-T/2) ] , 
a[pI(T/2)-pI(-T/2)]  +  6fpII(T/2)-pII(-T/2) ]  +  V(T/2)AT  *  -[pIII(T/2)-pIII( 
a[pI(T/2)-pI(-T/2)]  +  B[pII(T/2)-pII(-T/2)]  +  V(T/2)AT  =  -[pIII(T/2)-pIlI( 

At  this  stage,  we  should  notice  that,  along  a  symmetric  orbit,  the  normal 
displacement  n***  1j*  an  even  function  of  the  time  while  the  resulting 

tangential  displace! sent  p***  is  an  odd  function  of  the  time.  Therefore, 


(84) 
-T/2)  ] , 

-T/2)  ] . 
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the  predictor  formulae  (84)  become 

8nn(T/2)  =  0, 

anI(T/2)  =  -nIH(T/2), 

ap1  (T/2)  +  j  V(T/2)AT  =  -Pni(T/2), 

6pI(T/2)  *  0. 

As  we  have  seen  above,  if  we  exclude  the  critical  case  when  Tr(T)  *  2, 
we  must  take  8=0.  In  which  case,  the  second  of  equations  (85)  yields 
a  unique  determination  for  a  while  the  third  equation  provides  the 
tangential  coefficient  AT  =  2dTg/dC.  Thus  in  reference  to  a  symmetric 
c<  ..iodic  orbit, 


dc  vv.;  * 


"  2a  sin  $o’ 


dC  yo(co)  °* 


dc  xo(co)  °* 


dc  yo((V 


1  s*n  4>q  ^ 

2  “V~  +  2(2  Wa  Sln  ♦o' 


(86) 


We  should  observe  that,  from  a  symmetric  periodic  orbit  which  is  not 
singular  (Tr(T)  j  2),  the  tangential  predictor  generates  only  symmetric 
orbits.  In  other  words,  in  the  neighborhood  of  a  symmetric  periodic  orbit, 
a  natural  family  consists  only  of  symmetric  orbits,  and  it  can  be 
continued  by  a  sequence  of  nonsymmetric  ones  only  if  it:  goes  through  a 
symmetric  orbir.  whose  characteristic  exponent  is  zero. 
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8.  ISOENERGETIC  VARIATIONS  WITH  RESPECT  TO  A  PARAMETER 

We  consider  here  conservative  dynamical  systems  in  which  the 
characteristic  functions  A  and  W  depend  on  a  parameter  e.  Let  us 
assume  that  the  right-hand  members  of  the  equations  (1)  verify  conditions 
sufficient  for  the  solutions  of  (1)  to  be  unique  and  continuously 
differentiable  in  all  their  arguments  for  e  in  a  certain  neighborhood 
of  £q.  Suppose  that  x(t;e),  y(t;e)  is  a  particular  solution  of  (1) 
for  every  e  in  that  neighborhood  of  and  reduces  at  e^  to  the 

solution 

x(t),  y(t).  (87) 

Then  the  partial  derivatives 

u(t)  “  a7  x(t;e0),  v(t)  =  ~  y(t;eQ) 


are  solutions  of  the  variational  equations 

ii  *  2Av  +  (W  +2A  y)u  +  (W  +2A  y)v  +  (W  +2A  y)  , 
xx  xJ  xy  y  xe  ’ 

(88) 

v  =  -2Au  +  (W  -2A  x)u  +  (W  -2A  x)v  +  (W  -2A  x)  , 
xy  x  yy  y  ye  e  ’ 

to  which  belongs  the  variational  integral 


y=W  +Wu+Wv-xu-  yv,  (89) 

ex  y 

If  a  displacement  of  the  solution  (87)  caused  by  a  variation  Ae 
of  the  parameter  e  from  its  value  e^  for  the  solution,  i.e.  a 
solution  of  the  system  (88),  is  such  that  the  integration  constant 
defined  by  (89)  vanishes,  then  it  is  called  an  isoenergetio  variation 
of  the  solution  with  respect  to  the  parameter  e. 
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Functions  p(t)  and  n(t)  are  called  respectively  tangential 
and  normal  displacements  of  the  solution  (87)  with  respect  to  the 
parameter  e  if  the  equations  (88)  possess  a  solution  u(t),v(t) 
by  means  of  which  p(t)  and  n(t)  are  representable  in  the  form 
(15).  In  particular,  if  they  belong  to  a  displacement  u(t),v(t) 
for  which  the  integration  constant  y  in  (89)  vanishes,  then  they 
are  called  isoenergetic. 

Proceeding  exactly  as  in  paragraph  2,  one  is  able  to  prove  the 
following  proposition: 

THEOREM.  A  numerical  function  n  *  n(t)  is  an  isoenergetic 
normal  displacement  of  (87)  caused  by  a  variation  of  the  parameter  e 
if  and  only  if  it  satisfies  the  linear  differential  equation 

ft  +  On  -  -2W  ^  -2A  V  +  W  cos  <J>  -  W  sin  (90) 

e  V  e  ye  r  xe 


In  which  case >  a  numerical  function  p  -  p(t)  is  the  corresponding 
isoenergetic  tangential  displacement  of  (87)  if  and  only  if  it  is  given 
by  the  quadrature 


JL  (£) 

dt  V 


2n 


•  u 

A±i  +  _£ 

V  2 
V 


(91) 


The  homogeneous  linear  differential  equation  associated  to  (90)  being 
Hill's  equation,  the  isoenergetic  normal  displacements  of  (87)  with  respect 
to  the  parameter  are  evidently  of  the  general  form 

I  ,  e  II  ,  e 
n  ■  an  +  Bn  +  n  , 
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where  a  and  6-  are  constants  to  be  determined  by  the  initial  conditions, 
and  ne  is  the  particular  solution  of  the  nonhomogeneous  Hill's  equation 
(90)  satisfying  the  initial  conditioi  i 

n£  (0)  =  ne(0)  =  0.  (92) 

9.  ISOENERGETIC  CHANGE  OF  THE  TIME  VARIABLE 

The  practical  importance  of  the  theorem  stated  in  the  preceding 
paragraph  arises  from  its  applications  to  Celestial  Mechanics.  In  a 
number  of  problems  there  occur  fixed  singularities  which  are  not 
essential,  such  as  the  binary  collisions  in  the  Restricted  Problem. 

It  is  usual  to  remove  them  in  two  steps,  by  a  conformal  transformation 
of  the  coordinates  followed  by  a  change  of  the  time  variable  in  an 
isoenergetic  way.  Consequently,  in  the  transformed  system,  the 
original  constant  of  energy  turns  out  to  have  become  a  parameter,  and 
the  variations  with  respect  to  C  which  form  the  basic  ingredient  of 
our  tangential  predictor  are  deprived  of  any  meaning  for  the  original 
problem  unless  they  are  isoenergetic  in  reference  to  the  new  problem. 

To  convince  ourselves  how  essential  this  restriction  is,  let  us  recall 
what  is  meant  by  an  isoenergetic  change  of  the  time  in  Celestial 
Dynamics. 

Consider  a  conservative  dynamical  system  described  by  the  Hamiltonian 
function  Jo  in  the  phase  space  (q,p).  The  equations  of  motion 

q=£p,  P  -  -\>q  (93) 


admit  the  integral 
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$»(q,p)  =  h  (94) 

Denote  by  E(h)  the  manifold  defined  implicitly  by  (94)  in  the  phase 
space  (q,p).  If  only  solutions  lying  on  the  integral  manifold  E(h) 
for  a  fixed  value  of  the  integration  constant  h  are  considered,  one 
can  produce  in  a  direct  manner  a  rule  for  the  introduction  of  a  new 
time  variable. 

Let  G(q,p)  be  a  nonzero  numerical  function  defined  in  the  phase 
space  of  the  system.  Along  any  solution  (q(t),p(t))  of  the  equations 
of  motion  (93) ,  consider  the  new  time  variable  s  defined  by  the 
curvilinear  integral 

,  x  m  du _ 

S  J  G(q(u),p(u)) 

Define  a  conservative  Hamiltonian  function 

tf(q,p;h)  =  (§-h)G. 

Then  denoting  by  a  prime  differentiation  with  respect  to  s,  write  the 
canonical  equations  generated  by  X, 

q'  p'  =  -Xq. 

Denote  by  F(h)  the  manifold  defined  implicitly  by  the  invariant 
relation 

#(q»p;h)  =  0  (98) 


(96) 


(97) 


in  the  phase  space. 
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As  a  result  of  a  proposition  established  by  Painlev£,  only  the 
solutions  of  the  equations  (93)  in  the  integral  manifold  E(h) 
correspond  to  solutions  of  the  equations  (97),  and  these  lie  in  the 
manifold  F(h).  Conversely,  only  the  solutions  of  the  equations  (97) 
in  the  manifold  F(h)  correspond  to  solutions  of  the  equations  (93), 
and  these  lie  in  the  manifold  E(h) . 

This  duality  between  the  solutions  in  E(h)  and  those  in  F(h) 
is  to  be  extended  to  the  displacements  of  the  orbits. 

Of  a  displacement  (6q,6p)  of  an  orbit  (q(t),p(t))  solution  of 
the  equations  (93),  we  say  that  it  is  -isoenergetic  if  the  variation 
it  causes  to  the  Hamiltonian  function  is  zero.  In  the  same  way,  we  call 
^f-isoenergetic  a  displacement  of  an  orbit  (q,(s),p(s))  solution  of  the 
equations  (97)  such  that 


&X  -  0. 


The  orbit  obtained  by  displacing  an  orbit  lying  in  the  manifold  F(h) 
corresponds  by  Painlev^'s  duality  to  an  orbit  lying  in  E(h)  if  and  only 
if  the  displacement  is  ^(-isoenergetic.  In  which  case,  since 

6X  -  G  5$  +  (§  -h)  6  G  ■  G  S.'O 

in  view  of  the  fact  that  £  =  h  along  the  orbit,  to  a  displacement  which 
is  JK-isoenergetic  corresponds  by  duality  a  displacement  which  is 
JO-isoenergetic .  Conversely,  the  orbit  obtained  by  displacing  a  solution 
of  the  equations  (93)  lying  in  the  manifold  E(h)  corresponds  by  Painlev^'s 
duality  to  a  solution  of  (97)  lying  in  F(h)  if  and  only  if  the 
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displacement  is  ^-isoenergetic;  accordingly  to  a  displacement  which  is 
£>-isoenergetic  corresponds  by  duality  a  displacement  which  is 
JK-isoenergetic . 

In  consequence,  when  the  numerical  continuation  of  a  natural 
family  of  periodic  orbits  for  the  system  described  by  necessitates 
an  isoenergetic  change  of  the  time,  the  use  of  isoenergetic  displacements 
is  no  longer  optional,  but  it  is  a  necessity. 

In  the  phase  space  defined  by  the  Hamiltonian  X,  the  corrector 
part  of  our  procedure  goes  without  any  modification  at  all,  since  it  is 
based  on  isoenergetic  variations.  As  for  the  tangential  predictor,  we 
have  to  remark  that,  in  the  phase  space  of  X,  the  energy  constant  h 
is  to  be  treated  as  a  parameter.  Thus  the  rates  of  variation  n111  and 
p11*  which  we  can  compute  in  the  phase  space  $>  from  a  Hill  equation 
of  the  type  (16)  should  be  computed  in  the  phase  space  X  in  a 
^(-isoenergetic  way  from  a  Hill  equation  of  the  type  (90). 


v 
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CONCLUSIONS 

The  method  which  we  propose  here  for  continuing  numerically 
natural  families  of  periodic  orbits  is  valid  for  any  conservative 
dynamical  system.  But  to  date  we  have  applied  it  only  to  the 
planar  Restricted  Problem  of  Three  Bodies,  whether  in  barycentric 
synodical  Cartesian  coordinates  or  in  regularizing  parabolic 
coordinates.  The  program  has  generated  with  the  same  ease  both 
symmetric  and  nonsymmetric  orbits.  The  continuation  of  a  branch 
through  and  beyond  an  ejection  orbit  has  not  caused  any  hardship; 
the  calculation  of  characteristic  exponents  for  orbits  of  close 
approach  or  for  periodic  collision  courses  did  not  either  raise  any 
difficulty.  Results  in  these  directions  will  be  reported  soon. 

The  algorithm  involves  only  the  first  order  variations. 
Consequently,  it  is  unable  to  treat  the  singular  cases  of  periodic 
orbits  whose  characteristic  exponent  is  equal  to  zero.  These  appear 
as  singularities  on  the  analytic  manifold.  But  our  procedures  are 
open  to  the  introduction  of  variations  of  a  higher  order,  and  we  are 
presently  engaged  in  exploring  the  possibilities  offered  by  the  second 
variations.  This  way  we  hope  to  clarify  some  of  these  singular  orbits, 
at  least  those  at  which  the  natural  families  present  a  bifurcation  of 


order  2. 
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